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Abstract 

This  paper  is  an  overview  of  the  main  ideas  of  the  Generalized  Finite 
Element  Method  (GFEM).  We  present  the  basic  results,  experiences  with, 
and  potentials  of  this  method.  The  GFEM  is  a  generalization  of  the 
classical  Finite  Element  Method  —  in  its  h,p,  and  h-p  versions  —  as  well 
as  of  the  various  forms  of  meshless  methods  used  in  engineering. 

AMS(MOS)  subject  classifications.  65N30,  65N15,  41A10,  42A10, 
41A30 


1  Introduction 

A  numerical  method  to  approximate  the  solution  of  a  boundary  value  problem 
(BVP)  for  partial  differential  equations  (PDE)  has  two  major  components: 

(a)  The  selection  of  a  family  {ojj}jLi  of  small  sets  that  form  a  cover  of  the  do¬ 
main  of  the  BVP,  and,  for  each  j,  a  finite  dimensional  local  approximation 
space  Vj  of  functions  with  the  property  that  functions  in  Vj  can  accurately 
approximate  the  solution  of  the  BVP  on  u>j,  i.e.,  locally.  The  approximate 
solution  of  the  BVP  is  then  sought  from  the  space  S  of  global  functions, 
obtained  by  “pasting  together”  the  functions  in  Vj,  in  such  a  way  that 
good  local  approximability  of  the  Vj s  ensure  good  global  approximability 
of  S.  The  functions  in  S  are  often  of  the  form  'Yhj(t>jvji  with  Vj  £  Vj, 
and  where  {<fij}  is  a  partition  of  unity  with  respect  to  {u>j}.  We  note  that 
each  Vj  £  Vj  can  be  viewed  as  a  vector  of  real  numbers  (the  coefficients 
with  respect  to  some  basis  for  Vj).  Consequently,  a  functions  in  the  space 
S,  which  has  the  form  y~b  4>jVj,  can  also  be  viewed  as  a  vector  c  of  real 
numbers. 
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(b)  A  discretization  principle  that  selects  an  approximate  solution  of  the  BVP 
from  the  space  S;  in  other  words,  the  discretization  principle  associates  a 
specific  vector  c,  i.e.,  a  specific  element  of  S,  to  the  exact  solution  of  the 
BVP.  This  element  of  S  is  then  viewed  as  an  approximate  solution  of  the 
BVP. 

Given  the  local  spaces  Vj,  and  the  derived  global  space  S,  a  discretization  prin¬ 
ciple  determines  the  approximate  solution  in  S  by  approximating  the  partial 
differential  operator,  and  thereby  reduces  the  BVP  to  a  system  of  linear  or 
non-linear  algebraic  equations  for  the  vector  c.  When  the  system  is  linear,  the 
associated  matrix  is  often  sparse.  The  accuracy  of  the  approximate  solution 
depends  on  the  stability  of  the  discretized  partial  differential  operator  and  on 
the  approximation  properties  of  the  space  S,  which  in  turn  depends  on  the 
approximation  properties  of  the  spaces  Vj. 

We  first  discuss  briefly  the  choice  of  the  space  S,  as  indicated  in  (a),  for  dif¬ 
ferent  numerical  methods.  In  a  large  family  of  methods,  classical  interpolation 
theory  provides  guidance  in  the  choice  of  the  spaces  Vj,  and  thus  S.  Specifically, 
let  {xj}  be  a  set  of  given  distinct  points,  called  nodes ,  in  the  domain  of  defini¬ 
tion  of  the  BVP,  and  suppose  that  g  is  a  function  whose  values  gj  at  nodes  Xj , 
i.e.,  gj  =  g(xj),  are  given.  Then  the  space  S  is  such  that  there  exists  a  unique 
interpolating  function  f  £  S  such  that  f(xj)  =  gj.  The  approximation  property 
of  the  space  S  is  related  to  the  interpolation  error,  i.e.,  g  —  f,  and  this  error 
depends  on  the  distribution  of  the  nodes  {ay},  which  could  be  regular  or  irregu¬ 
lar  (scattered  nodes),  and  on  the  bounds  of  higher  derivatives  of  the  function  g. 
The  space  of  polynomials,  piecewise  polynomials,  and  the  combination  of  radial 
basis  functions  are  examples  of  the  space  S  with  this  interpolation  property.  We 
mention  that  the  uniqueness  of  the  interpolating  function  f  £  S,  with  respect 
to  the  given  data  {gj},  depends  strongly  on  the  distribution  of  nodes,  as  well  as 
on  the  space  S  (and  thus  on  spaces  Vj).  For  a  given  distribution  of  the  nodes, 
the  space  S  may  not  have  unique  solvability  of  the  interpolation  problem.  To 
resolve  this  problem  in  certain  situations,  various  stabilization  techniques  have 
been  reported  in  the  literature;  e.g.,  see  [10]  in  the  context  of  thin-plate  spline 
radial  functions.  We  also  refer  to  [11]  for  a  detailed  discussion  on  radial  basis 
functions.  The  interpolation  problem  for  the  space  of  polynomials  or  piecewise 
polynomials,  and  its  sensitivity  on  the  distribution  of  nodes,  is  well  studied  in 
the  literature. 

But  there  are  other  methods,  e.g.,  certain  meshless  methods,  where  the 
choice  of  Vj,  and  thus  S,  is  not  dictated  by  the  idea  of  interpolation.  In  these 
methods,  the  local  spaces  Vj  are  constructed  from  particle  shape  functions,  e.g., 
RKP  shape  functions,  and  the  elements  of  the  space  S  are  of  the  form  Yhjvj, 
where  v3  £  Vj  {i.e.,  <f)j  =  1  in  Y}j(f>jvj)-  The  approximability  of  the  spaces 
Vj  and  S,  is  ensured  by  so  called  “reproducing  property”  of  the  particle  shape 
functions.  For  a  detailed  discussion  of  these  spaces,  we  refer  to  [4,  26]. 

We  will  now  briefly  discuss  the  discretization  principle  indicated  in  (b).  Dif¬ 
ferent  discretization  principles,  together  with  given  global  approximating  spaces 
S,  give  rise  to  different  methods  for  the  approximation  of  the  solution  of  a  BVP; 
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e.g.,  finite  difference  methods  (FDM),  finite  volume  methods  (FVM),  colloca¬ 
tion  methods,  and  methods  based  on  weighted  residuals.  We  note  that  the  FDM 
and  collocation  methods  can  be  viewed  as  obtained  from  the  discretization  by 
the  Petrov-Galerkin  method  (in  the  most  general  setting)  with  Dirac  functions 
used  as  test  functions.  Establishing  stability  and  obtaining  error  estimates  for 
these  methods  is  subtle  and  difficult,  even  when  the  spaces  Vj,  and  consequently 
S,  have  good  approximation  properties.  For  example,  though  the  convergence 
analysis  of  FDM  with  regularly  distributed  nodes  is  well-established,  not  much 
is  known  when  the  nodes  are  irregularly  distributed  ([14]).  The  convergence  of 
the  collocation  method  using  radial  functions  was  analyzed  in  [18].  For  a  survey 
of  application  of  these  methods,  we  refer  to  [24], 

A  variant  of  the  collocation  method,  obtained  from  the  discretization  by  the 
Petrov-Galerkin  method  using  test  functions  with  small  supports  (instead  of 
Dirac  functions  as  mentioned  in  the  last  paragraph),  have  also  been  reported  in 
the  engineering  literature,  but  without  rigorous  mathematical  analysis  ([1,  25, 
43,  32]).  These  methods,  which  are  also  used  to  approximate  solutions  of  non¬ 
linear  equations,  CFD,  and  other  engineering  problems,  often  lack  robustness. 
Various  ad-hoc  stabilization  techniques  are  used  in  the  implementation  of  these 
methods,  without  rigorous  mathematical  examination. 

There  is  yet  another  class  of  methods  that  is  based  on  a  discretization  prin¬ 
ciple  where  the  trial  and  test  functions  belong  to  the  same  Hilbert  space,  say 
the  Sobolev  space  Ff1(D)  (for  second  order  elliptic  BVP).  This  principle  is  re¬ 
ferred  to  as  the  Galerkin  method  or  Bubnov-Galerkin  method  ([30]).  Typical 
representatives  of  this  class  are  Finite  Element  Method  (FEM)  -  with  its  h,  p, 
and  h  —  p  versions  and  mixed  FEM.  In  these  methods,  the  functions  in  the  local 
spaces  Vj  have  to  be  “pasted  together”  so  that  the  space  £  is  a  subspace  of 
U1(D).  This  is  achieved  by  considering  Vj’ s  consisting  of  piecewise  polynomials 
(or  pull-back  polynomials)  of  special  form,  defined  with  respect  to  an  appropri¬ 
ate  mesh.  Certain  classes  of  meshless  methods,  e.g.,  RKP  method,  fall  in  this 
category.  In  these  meshless  methods,  the  spaces  Vj  are  subspaces  of  the  energy 
space,  and  consequently  the  elements  of  S,  which  are  linear  combinations  of 
elements  in  Vj  (mentioned  before),  are  automatically  in  the  energy  space. 

In  this  paper,  we  present  the  main  ideas  of  Generalized  Finite  Element 
Method  (GFEM),  which  is  a  Galerkin  (or  Bubnov-Galerkin)  method.  The  lo¬ 
cal  spaces  Vj  consists  of  functions,  not  necessarily  polynomials,  that  reflect  the 
available  information  on  the  unknown  solution  and  thus  ensure  good  local  ap¬ 
proximation.  Then  a  partition  of  unity  {<Pj}  is  used  to  “paste”  these  spaces 
together  to  form  S,  which  is  a  subspace  of  the  energy  space  and  has  good  global 
approximability.  The  GFEM  has  been  extensively  discussed  in  a  series  of  papers 
([16,  37,  38,  40,  39]),  where  its  effectiveness  was  shown  when  applied  to  prob¬ 
lems  with  domains  having  complicated  boundaries,  problems  with  micro-scales, 
and  problems  with  boundary  layers.  We  will  present  the  theoretical  basis  of  the 
GFEM,  proving  major  results.  In  addition,  we  will  discuss  various  procedures 
for  the  selection  of  local  approximating  functions  and  comment  on  certain  issues 
in  implementation. 
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The  partition  of  unity  approach  was  first  used  in  [5]  to  obtain  an  accurate 
approximation  to  the  solution  (which  is  non-smooth)  of  BVP  for  PDEs  with 
rough  coefficients;  the  method  in  [5]  was  referred  to  as  the  Special  Finite  Element 
Method.  The  importance  of  such  an  approach  was  seen  in  [7],  which  showed  that 
standard  FE  approximations  converge  arbitrarily  slowly  when  approximating 
solutions  to  problems  with  rough  coefficients.  Based  on  the  ideas  in  [5],  the 
GFEM  was  elaborated  on  in  [6,  27,  28],  where  it  was  referred  to  as  the  Partition 
of  Unity  Method  (PUM).  Later  in  [37,  38],  the  method  was  referred  to  as  GFEM, 
since  the  classical  FEM  is  a  special  case  of  this  method.  Currently,  the  partition 
of  unity  approach  is  used  in  various  directions  under  various  names  —  Method 
of  Clouds,  XFEM  (extended  FEM),  and  Method  of  Spheres  ([15,  41,  36,  13]). 
These  methods  differ  primarily  in  the  form  of  partition  of  unity  functions  used 
and  in  the  use  of  different  local  spaces. 

2  The  Galerkin  Method 

Suppose  we  are  interested  in  solving  the  stationary  heat  conduction  problem  on 
the  domain  fi  C  R2  with  piecewise  smooth  boundary  r  =  Ti  ur2.  Specifically, 
we  consider  the  problem 

!— div  (a(x,  y)  grad  u )  =  /,  for  (a:,  y)  €  fi 

«  =  0onFi  (2.1) 

a§F  =  5  on  r2. 

Here  /  =  f(x,y)  is  the  heat  gain  from  internal  sources  per  unit  volume,  a  = 
a(x,  y)  is  the  coefficient  of  thermal  conductivity,  g  =  g{x,  y)  is  the  heat  flow 
per  unit  length  across  r2.  We  consider  /,  a,  and  g  to  be  given,  we  specify  the 
temperature  to  be  0  on  F-| ,  and  specify  the  heat  flow  per  unit  length  across  T2  to 
be  g ,  and  seek  the  steady  state  temperature  u  =  u(x,  y)  throughout  the  domain 
fL  The  function  a(x,y)  could  be  rough,  i.e.,  fail  to  have  continuous  derivatives, 
but  is  assumed  to  satisfy 

0  <  a  <  a{ x,  y)  <  j3  <  oo. 

As  usual,  we  give  our  problem  a  weak,  or  variational,  formulation.  Let 

£(fl)  =  £  =  {v  :  |M||(n)  <  oo | ,  (2.2) 

where 

IMli(fi)  =  IMI  e  =  J^a(x>y)  (jj^j  +  (|^)  dxdy .  (2.3) 

We  note  that  under  mild  restrictions  on  T,  ||u||f  (n)  <  00  implies 

IMl!2(n)  =  IMl!2  =  [  a\u\2dxdy  <  oo,  (2.4) 

“  a  Jn 
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i.e.,  v  G  L*(fi).  We  then  let 


£r1  (f2)  =  £n  =  {v  :  v  G  £(fl),  u  =  0  on  I\}  ,  (2-5) 


where  the  Dirichlet  boundary  condition  is  imposed  in  the  sense  of  trace.  If 
Id  =  0,  then  £ri(f^)  =  £(f2).  The  space  £r,  is  the  energy  space  for  our  problem 
and  ||f  ||f  is  the  energy  norm  of  v.  (Strictly  speaking,  ||f  ||f(o)  is  not  a  norm  on 
£(fi);  it  is,  however,  a  norm  up  to  rigid  body  motions,  which  in  this  situation 
are  the  constants.) 

On  £r  ,  x  £r1  define 


B(u,  v ) 


du  dv 
dx  dx 


du  dv 
dy  dy 


dx  dy 


and 


F(v)  =  /  fvdxdy+  /  gvds, 

J n  '  J r2 

where  we  assume  /  G  L2(£l)  and  g  G  L2(T2).  Then  the  weak  formulation  reads, 


J  Find  u  G  £r1  satisfying 
\  B[u,v)  =  F{v)  for  all  v  G  £rv 


Remark  2.1.  If  T2  =  0,  then  (2.1)  is  Dirichlet  Problem.  If  the  lengths  of 
both  Ti  and  T2  are  positive,  then  (2.1)  is  a  mixed  Dirichlet-Neumann  Problem. 
In  either  of  these  cases,  Problem  (2.1)  ((2.6))  is  uniquely  solvable.  If  P i  =  0, 
then  (2.1)  is  a  Neumann  Problem.  In  this  case  (2.1)  ((2.6))  will  be  solvable 
provided  Jn  fdxdy  +  fr  gds  =  0.  To  ensure  uniqueness,  one  needs  an  auxiliary 
condition:  say  fn  udxdy  =  0. 

We  next  consider  the  approximation  of  the  solution  u  of  (2.1)  ((2.6))  by  the 
Galerkin  Method  (Bubnov-Galerkin  method).  Toward  this  end  we  suppose  we 
have  a  finite  dimensional  space  S  C  £r, ,  and  consider  the  problem 

J  Find  us  G  S  satisfying  ,  , 

\  B(us,v)  =  F{v)  for  all  v  G  S. 

This  problem,  like  Problem  (2.6),  has  a  unique  solution,  and  is  equivalent  to  a 
system  of  linear  algebraic  equations.  Specifically,  ii  <fii, ... ,  4>m  spans  the  space 
S  and  we  write  us  =  Sjli  Problem  (2.7)  becomes 

n 

B{<j>i,  <t>j)c:j  =  F((j>i),  i  =  1, . . . ,  m.  (2.8) 

j= i 

If  is  a  basis  for  S,  then  the  linear  system  (2.8)  is  nonsingular  and  is 

uniquely  solvable.  If  {(f>j}1JL1  is  not  a  basis,  i.e.,  it  fails  to  be  linearly  indepen¬ 
dent,  the  system  (2.8)  is  solvable  since  (2.7)  is  solvable,  but  solutions  of  (2.8)  are 
not  unique.  (The  family  {<pj}JL1  is  said  to  span  S  if  any  v  G  S  can  be  written 
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as  v  =  J2jLi  Cj(f>j  for  some  coefficients  Cj ;  it  is  said  to  be  a  basis  if,  in  addition, 
it  is  linearly  independent,  i.e.,  J2jLi  Cjrf’j  =  0  implies  Cj  =  0,  j  =  1, . . . ,  m).  We 
note,  however,  that  if  {cj1)}JL1  and  {c^2)}fL1  are  solutions  of  (2.8),  then 

m  m 

i-l  j=i 

Whenever  we  have  a  spanning  set  <f> i, . . . ,  4>mi  we  refer  to  the  functions  (j)j  as 
shape  functions.  If  the  shape  function  have  small  supports,  the  matrix  of  the 
system  (2.8)  is  sparse.  The  Finite  Element  Method  (FEM)  is  of  this  type,  with 
piecewise  polynomial  shape  functions  defined  on  a  mesh. 

We  will  measure  the  accuracy  of  us  in  the  energy  norm.  Letting  es  =  u  —  us 
be  the  error ,  and  consider  the  energy  norm  of  the  error. 

\\es\\e  =  (B(es,es))1/2.  (2.9) 

The  main  feature  of  the  Galerkin  Method  is  that 

||«-  usIU  =  lies’ll  <  ||u-£||f,  for  any  f  £  S.  (2.10) 

We  thus  need  to  construct  S  so  that 

Sc£ Fl(fi)  (2.11) 

and  so  that 

there  exists  (  =  €  S  so  that  ||u  —  £“||f:(fi)  is  small  .  (2-12) 

Of  course,  it  is  also  important  that  the  approximating  space  S  lead  to  a  reason¬ 
ably  solvable  linear  system  (2.8).  Constructing  S  so  that  (2.11)  and  (2.12)  are 
satisfied  are  our  major  goals. 

In  many  important  problems  the  character  (smoothness)  of  the  solution 
changes  from  one  part  of  the  domain  to  another,  so  it  is  natural  to  attempt 
to  approximate  u  separately  on  these  parts  of  O.  There  is  often  a  natural  di¬ 
vision  of  f l  into  subdomains,  uij,  so  that,  for  each  j,  we  can  find  a  function  Eff 
that  approximates  u  well  on  cur  More  precisely,  we  have  open  sets  u> i, . . . ,  wjv, 
called  patches  satisfying 


N 

u>j  C  O  and  Q  =  luj  (they  form  an  open  cover  of  f 1),  (2-13) 

i= i 

and  function  G  £(tOj)  satisfying 

\\u  -  is  small,  (2.14) 

where  £(uij)  and  ||u  —  \\e(u>j)  are  defined  by  (2.2)  and  (2.3),  with  replaced 

by  uij .  We  will  speak  of  {ujj}  as  a  partition  of  O.  We  then  need  to  “paste”  these 
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approximating  functions  together  to  obtain  a  function  G  S  satisfying  (2.12). 
These  two  aspects  of  our  development  —  the  existence  of  local  approximations 
and  the  process  of  pasting  them  together  -  are  largely  independent. 

For  each  j  we  wish  first  to  construct  £“  on  u>j  so  that  (2.14)  is  satisfied. 
Then  we  wish  to  construct  using  the  £“  —  pasting  them  together  —  so 

that 

N  N 

Ki  E  ii«  -  ^iiifo,,)  <  ii«  -  ni(o)  <  k2  e  ii«  -  (2-!5) 

t=i  1=1 

where  K\,K2  are  independent  of  u  and  the  number  of  patches  (N),  but  do 
depend  on  the  form  (character)  of  the  patches.  Our  main  focus  in  Section  3  will 
be  to  prove  the  upper  bound  in  (2.15).  The  lower  bound  will  be  true  in  some 
situations,  but  not  in  others. 

These  issues  will  be  discussed  in  detail  in  the  next  section.  We  end  this 
section  by  noting  that  to  find  a  suitable  £“  and  to  show  that  ||w  —  £ j\\s(wj ) 
is  small,  we  need  to  use  the  available  information  on  the  (unknown)  solution. 
For  example,  if  a(x,y),  f(x,y),  g(x,y)  are  smooth  functions  and  Ti  =  T  is  also 
smooth,  then  u(x,  y)  will  be  a  smooth  function.  From  standard  polynomial 
approximation  theory  we  thus  know  that  there  is  a  quadratic  polynomial 
that  approximates  u  well  on  uif 

l«  c  •'/';/'/ > 

where  Kj  is  a  bound  on  the  third  derivatives  of  u  on  ujj  (\D3u\  <  Kj)  and  hj  is 
the  diameter  of  u>j  and  Cj  depends  on  the  form  of  u)j. 

3  Local  and  Global  Approximation 

In  this  section  we  show  how  to  accomplish  the  goals  stated  in  Section  2  - 
namely  (2.11)  and  (2.12)  —  by  means  of  local  approximation  and  the  pasting 
process,  which  are  largely  separate.  As  indicated  in  Section  2,  let  {ujj}f=1  be 
open  sets  (patches)  satisfying 


N 

cdj  C  ff  and  O  =  |^J  tOj. 

3=  i 

We  assume  in  addition  that  any  ie!l  belongs  to  at  most  n  of  the  subdomains 
ojj.  Then  let  {<f>j}jL i  be  a  family  of  functions  defined  on  Cl,  having  piecewise 
continuous  first  derivatives,  and  satisfying  the  following  properties: 

<i>j(x,y)  =  o,  for  (x,y)  G  Cl  \  loj,  j  =  1, . . . ,  N-  (3.1) 

N 

E  y) =  for  (x’V) e  (3-2) 

3= 1 
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(3.3) 


max  \<f>j{x,  y)\  <  C\,  j  =  1, . . . ,  N\ 
(x,y)en 


max  \X7(j)j(x,  y)\  < 

(x,y)eQ. 


C2 

diam  (uj) 


(3.4) 


where  0  <  C\ ,  C'2  <  00.  Here  diam  (uj)  denotes  the  diameter  of  ojj .  Property 
(3.2)  states  that  is  a  partition  of  unity  on  fi. 

As  an  example,  consider  the  classical  FEM  with  triangular  elements  satis¬ 
fying  the  minimal  angle  condition,  with  nodal  points  Aj.  Let  Uj  be  the  patch 
or  finite  element  star  associated  with  the  node  Aj,  i.e.,  the  union  of  triangles 
with  Aj  as  one  of  their  vertices.  It  is  easy  to  see  that  the  family  ojj  creates  a 
partition  of  f 1.  Further,  let  <f>j  be  the  piecewise  linear  functions  with 


4>i  (Aj ) 


1,  if  i  =  j 
0,  if  i^j. 


Then  it  is  easily  seen  that  the  family  {<f>j}  satisfies  (3.1)-(3.4)  with  C\  =  1  and 
C'2  depending  on  the  minimal  angle  condition. 

We  next  mention  another  example.  Let 


ft  =  {(#,  y)  :  0  <  x  <  1, 0  <  y  <  1} 

and  let  A^  =  Ai  j  =  ( ih,jh),h  =  —,i,j  =  0,1,.  ,.,m.  Let  lu £  be  the  intersec¬ 
tion  of  H  and  the  open  disk  centered  at  A/.  with  radius  Rh ,  where  R  is  such 
that  {u>k}  is  a  cover  of  O.  Letting  4>(r),  0  <  r  <  00,  be  a  function  with  bounded 
first  derivative  and  with  <j>(r)  >  0,  for  0  <  r  <  R,  and  4>(r)  =  0  for  r  >  R ,  define 


^k\x,y)  =  <t> 


The  family 
we  define 


satisfies  (3.1),  (3.3),  and  (3.4),  but  not,  in  general,  (3.2).  If 


k(x,y ) 


y) 

E  i4>i(x,vY 


then  the  family  {</>%}  satisfies  all  the  conditions  (3.1)-(3.4).  To  prove  (3.3)  and 
(3.4)  for  this  family,  we  use  the  fact  that 


^2  & ih) (x’y)^T>°i  for  (x> V)  e 

1 


The  functions  in  the  family  {</)%}  are  called  Shepard  functions  ([35]). 

To  every  ujj  of  the  partition  {uij}  we  associate  an  rn  (j )  -  d  i  mensi  on  a!  space  of 
functions  defined  on  u>j : 

m(j) 

1)  —  {0  ■  ^  [  bji£ji,  bji  (E  ft ,  fj i  C  £(wj)  and  fji  —  0  on  ajj  C  I  | }  •  (3.5) 

i= 1 
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The  space  Vj  is  called  a  local  approximation  space.  Note  that  the  (essential) 
Dirichlet  boundary  condition  is  built  into  Vj.  Then  we  let 


gGFEM  =  =  4>jtj  ■  where  &  G  Vj  \ 

=  span  of  {rjji,  i  =  1, . . . ,  =  1,...,N},  (3.6) 

where 

Vji  =  (3-7) 

are  the  shape  functions  for  the  SGFEM .  The  space  SGFEM  is  called  the  Gen¬ 
eralized  Finite  Element  global  approximation  space. 

Theorem  3.1  We  have 

SGFEM  csTl(n).  (3.8) 

Proof.  Using  (3.1)  we  see  that  {<pjtiji)(x,y)  =  0  for  (x,y)  G  dutj  fl  ft.  Hence 
(fjf-jt  can  be  extended  as  zero  to  all  of  fl,  and  (pjfji,  so  extended,  will  be  in 
£(fi).  Furthermore,  since  fji  =  0  on  Uj  n  r1?  we  see  that  (pjZji |rx  =  0-  So,  for 
all  j  and  i,  G  £n(fi)>  and  hence  the  span  of  these  functions  is  in  £(fl). 

This  is  the  desired  result.  □ 

Remark  3.1.  Theorem  3.1  establishes  (2.11),  one  of  the  goals  discussed  in 
Section  2. 

The  Generalized  Finite  Element  Method  (GFEM)  is  now  defined  to  be  the 
Galerkin  Method  (2.7)  with 

g  —  gGFEM 

We  denote  the  approximate  solution  by  us  =  ugfem  ■  If  we  can  now  construct 
a  £“  G  SGFEM  so  that  (2.12)  is  satisfied,  then  from  (2.10)  we  know  that  ||u  — 
uGfem ||£  jg  sman  yye  now  turn  to  the  construction  of  such  a  fu. 

For  each  j,  we  assume  the  exact  solution  u  of  Problem  (2.1),  more  generally 
any  u  G  £rn  can  be  accurately  approximated  on  u>j  by  a  function  £“  G  Vj ; 
specifically  that 


(3.9) 

(3.10) 

(3.11) 
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We  see  that  the  local  approximation  is  ensured  by  the  appropriate  selection 
of  the  spaces  V);  and  the  pasting  together  is  handled  by  multiplication  by  the 
partition  of  unity  functions,  <j>j.  We  now  estimate  ||u— C]t||i2(n)  and  ||u— ||£(fi)- 


Theorem  3.2  Suppose  u  e£ri(U).  Then 


N 


1/2 


\u  ~  Chi(n)  <  k1/2Ci  (  4 O') 

ki=i 


(3.12) 


and 


m 


N 


1/2 


■-rikn)<(2  K)1/2 

1  :  aiann  (ujj)  .  n 

3  =  1  v  J  J  j  —  i 


+  Cf^2elU)  >  (3.13) 


where  C\  and  C2  are  as  in  (3.3)  and  (3-4),  respectively. 

Proof.  We  will  first  prove  (3.12).  Recalling  the  definition  of  fu  in  (3.11)  and 
using  the  fact  that  {(f>j}  is  a  partition  of  unity  on  Cl,  we  have 


\\u-f 


C\2dxdy 


(3.14) 


Using  the  fact  that  any  x  €  Cl  is  in  at  most  k  subdomains  ujj  we  see  that  the 
sum  ((j(u  ~  £3)  has  a^  most  k  terms  for  any  (x,  y)  €  Cl.  Hence,  using  the 

Schwartz  inequality,  we  have 


N  N 

l£fc(«-S,)la<"£lfc(«-tf)la- 

1=1  1=1 

Thus,  using  (3.3)  and  (3.9)  in  (3.14),  we  have 


2 

L2(f2) 


< 


< 


which  is  (3.12). 


(3.15) 


10 


Now  we  turn  to  the  proof  of  (3.13),  which  is  similar.  Proceeding  as  above, 
we  have 


I  u-£ 


v  II  2 


f  a|V  (u  —  £“)|2  dx  dy 

J  n 

N 

/  alV^faiu-^fdxdy 


l=i 

N 


a\  _  £l  )V<Al  +  <Pj^(u  -  £“)]| 2dxdy 


l=i 


<  2  ^  a  )V^j|  dx  dy  +  2  a  ^^cj)jV(u  —  dxdy 

..  N  ..  N 

<2  k  a  E  | («  -  C“) f  dxdy  +  2K  aJ2  |0lV(u  -  |2  dxdy. 

V  UJj  _ 1  j  Oil  _ 1 


<0  j=i 

Hence,  using  (3.9)  and  (3.10),  we  obtain 


^  j=i 


iV  2/0  JV 

-  nu  <  2«  ( cf  £  £l(2J)  +  c?  £  4c/)  ] , 

1  diam  (u;?  )  ' 

j— i  v  j=i 


(3.16) 


which  is  (3.13).  □ 

Since  £2(1)  is  usually  proportional  to  £i(j)/  diam  (cCj),  the  terms  in  (3.13) 
are  in  some  sense  balanced.  The  next  theorem  gives  sufficient  conditions  to 
ensure  this  balance. 


Theorem  3.3  Suppose  u  G  £rt  ■  Suppose  the  patches  {uij}  and  the  local  ap¬ 
proximation  spaces  {Vj}  satisfy  the  following  assumptions: 

(a)  For  all  j  for  which  uJj  D  Id  =  0,  Vj  contains  the  constant  functions,  and 


IMI Ll(wj)  <  Czdiam^j^vWs^.),  for  all  v  G  £(uj)  satisfying  /  avdxdy  =  0, 

J  UJj 

(3.17) 

i.e.,  for  all  v  with  weighted  a— average  over  u>j  equal  to  0; 

(b)  For  all  j  for  which  \u>j  nTi|  >0, 

\\v\\ Ll{uj)  <  C4diam(ajj)\\v\\£(u,.'),  for  all  v  G  £(u)j)  with  n^-nri  =  0.  (3.18) 

(Note  that  we  require  C3  and  C4  to  be  independent  of  j).  Then  there  exists 
G  Vj  so  that  the  corresponding  global  approximation, 

N 

r  =  EMl’  (3-19) 


1=1 
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satisfies 


1/2 


||«  -  ?“IU=(n)  <  C5  diam2(u)j)el(j)^  ,  (3.20) 

where  C5  =  y/RC\{C^  +  C2)1^2,  and 

N 

\\n-C\\£<C6(J24(j))1/2,  (3.21) 

j= 1 

where  C6  =  {2 K  {C\  +  C|(Cf  +  Cl))}1/2. 

Remark  3.2.  Estimates  (3.17)  and  (3.18)  are  Poincare  inequalities.  In  Re¬ 
marks  3.4  and  3.5,  we  give  simple  geometric  conditions  on  the  patches  uij  that 
imply  (3.17)  and  (3.18)  hold  uniformly  in  j.  Specifically,  we  bound  C3  and  C4 
in  term  of  simple  geometric  data. 

Proof.  Let  £“  satisfy  (3.9)  and  (3.10).  We  divide  the  index  set  A  = 
{1, . . . ,  N}  into  two  disjoint  sets: 

Air*  =  {j  ■■  1  <  j  <  N,uj  n r4  =  0} 

and 

Abd  =  {j  ■  1  <  j  <  N,  Uj  n  Ti  ^  0}. 

For  j  G  Alnt,  let  £“  =  ff  +  rj,  where  r,  is  a  constant  chosen  so  that  u  —  ff 
has  zero  a-average  on  uij.  By  assumption  (a),  G  Vr  Then,  using  (3.17)  with 
v  =  u  —  £“  and  noting  that  V(u  —  £“)  =  V(u  —  £“),  from  (3.10)  we  have 

\\u~  <  C2diam2K)  f  a\V(u-^)\2dxdy 

Jujj 

=  C|diam2(wj)  f  a\S7 (u  —  £fj)\2 dx  dy 

Jujj 

<  C2  diam2(wj)  e|(j).  (3.22) 

We  also  have 

ll'u^l“ll£(a;J)  =  [  a\V(u-gt)\2dxdy  <  e2(j).  (3.23) 

Jujj 

For  j  G  Abd,  let  =  £“.  Now  u|^nri  =  0,  and  we  know  that  fjl^nri  =  0. 
Thus,  using  (3.18),  with  v  =  u  —  £“,  and  (3.10),  we  have 

=  ||«— £“||l2(^)  <  C4  diam(wj)||u-^“||£(a;.)  <  C4diam(wj)e2(j)- 

(3.24) 

Also, 

\\u-ii\\eM<4U)-  (3-25) 
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Following  the  steps  leading  to  (3.15)  in  the  proof  of  Theorem  3.2,  and  using 
(3.22)  and  (3.24),  we  get 

iiu _  Iuiil2(n)  ^  k^i  iiM ~ iyiiiitui) 

jeA 

=  K(^i  |  X/  Hu  —  £j  Hu  _  IlijK) 

j&Abd 

<  ^(Cf  +  CSj^diamV^y),  (3.26) 

jeA 

which  is  (3.20)  with  C5  =  \fnC\  ( C'f  +  C2)1/2.  Similarly,  following  the  steps 
leading  to  (3.16)  in  the  proof  of  Theorem  3.2,  and  using  (3.22)-(3.25)  we  obtain 


llu  ~  f“lll  <  2n(Cl  +  C|((7|  +  C4))  ^  e|(i))  (3.27) 

jeA 


which  is  (3.21)  with  C6  =  y/2 k(C|  +  Cf  (Cf  +  Cl))1/2.  □ 

The  idea  of  GFEM,  in  particular  the  use  of  a  partition  of  unity  and  local 
shape  functions,  was  first  introduced  in  [5].  A  result  similar  to  Theorems  3.2 
and  3.3  was  proved  in  that  paper.  The  GFEM  was  further  developed  in  [6,  28]. 
Our  presentation  of  Theorems  3.2  and  3.3  closely  follows  [6,  28]. 

Remark  3.3.  Theorem  3.3  establishes  (2.12),  the  second  goal  discussed  in 
Section  2. 


Remark  3. 4-  Suppose  each  ujj  is  convex,  dj  =  diain(cCj),  and  ujj  contains  a 
ball  of  diameter  dj  >  ,  with  k-\  independent  of  j.  Then 


C3  ^  2ki 


(3.28) 


where  C3  is  the  Poincare  constant  in  (3.17).  This  follows  directly  from  Theorem 
8.1  in  the  Appendix  (Section  8). 

Remark  3.5.  Suppose  ujj  fl  Ti  is  an  arc.  Let  SA-nri^)  be  the  sector  sub¬ 
tending  this  arc,  and  let  7^.nri  be  the  angle  of  this  sector.  Suppose  each  Uj  is 
convex,  dj  =  diam  (oJj),  and  suppose  u>j  is  a  disk  of  diameter  dj  >  whose 
closure  lies  in  u>j.  Assume 

JujnrAx)  >  70,  for  all  x  €  u)j,j  =  1,  2, ...  ,1V. 


Then 


2*,+  |  t)™ 

'  J  7o 


(3.29) 


where  C4  is  the  Poincare  constant  in  (3.18).  This  follows  directly  from  Theorem 
8.2  in  the  Appendix  (Section  8). 
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Remark  3.6.  In  Theorems  3.1  and  3.2  we  have  imposed  only  minimal  condi¬ 
tions  on  the  patch  u>j.  In  Theorems  3.3  we  imposed  additional  conditions.  We 
note,  however,  that  the  conditions  on  the  u)j  can  be  considerably  relaxed.  The 
u)j  can,  in  particular,  be  multiply  connected.  The  condition  that  uJj  fl  Tj  is  an 
arc  can  be  relaxed;  in  particular,  it  can  be  a  disconnected  set  (see  Remark  8.1). 

We  return  now  to  the  GFEM.  Suppose  the  hypotheses  of  Theorem  3.3  are 
satisfied,  and  suppose  u  is  the  solution  of  (2.6).  It  follows  from  (2.10),  with 
£  =  £“,  and  (3.21)  that 

1  /2 

\\u  -  ugfem\\s(Q)  <  c\\u  -  <  c  fy^ejQ'))  ,  (3.30) 

which  is  the  main  error  estimate  for  the  GFEM.  It  will  be  useful  to  state  this 
estimate  in  the  following  alternate  form: 

l/2 

||«  -  ugfem llfyn)  <  C^inf^  (E  \\u  -  •  (3.31) 

We  can  write  uGFEM  as 


N  m(j) 

UGFEM  =  EE  CjiVji,  (3.32) 

l=i  »= 1 

where  c  =  {cjt}  is  the  solution  of  the  linear  system  (see  (2.8)) 

N  m(j) 

EE  Vji)cji  =  F(vik),  1  <  k  <  m(l),  1  <  l  <  N, 

l=i  »= 1 
or 

Ac  =  F,  (3.33) 

where  A  is  the  stiffness  matrix,  whose  elements  are 


A(l,k;j,i)  =  B(riik,r]ji)  =  /  ^Vik  •  Vrjjidx,  (3.34) 

J  ujj  n  u>i 

and  F  is  the  load  vector,  whose  components  are 

F(l;  k)  =  (  frjik  dx  +  f  grjik  ds.  (3.35) 

J  uji  JTi  Ouji 

The  GFEM  is  a  very  general  method.  We  show  in  the  next  section  that 
it  is  an  umbrella  covering  many  standard  FEMs,  hence  the  name  Generalized 
FEM.  Using  polynomial  functions  together  with  other  special  functions  we  get 
the  XFEM  (see  [36,  41]),  which  is  a  special  case  of  the  GFEM.  The  specific 
selections  of  (j>j  and  V3  lead  to  the  methods  referred  to  in  the  literature  by 
different  names. 
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Remark  3.6.  We  have  addressed  only  second  order  boundary  value  problems. 
In  an  analogous  way  the  GFEM  can  be  used  to  approximate  the  solutions  of 
2mth  order  boundary  value  problems,  where  the  bilinear  form  includes  deriva¬ 
tives  of  orders  up  to  m.  Instead  of  (3.4)  we  would  assume 


max 

(x,y)en 


\Da^j{x,  y)\  < 


C2 

(diam  C0j)m' 


where  a  =  (l,k),l,k  >  0 ,k+l  =  m.  In  addition,  we  have  to  assume  that  <j)j 
has  piecewise  continuous  derivatives  of  orders  up  to  m  on  Cl,  and  that  (j>3  and 
its  normal  derivatives  of  orders  up  to  m  are  0  on  uj3  n  Ti . 


4  Relation  Between  GFEM  and  Classical  FEM 

The  GFEM  is  based  on  the  generalization  of  the  idea  of  classical  FEMs.  We 
will  illustrate  this  by  showing  that  certain  classical  FEMs  can  be  cast  in  the 
framework  of  a  GFEM  by  appropriately  choosing  the  partition  of  unity  functions 
{4>j}  and  the  local  approximation  spaces  { V3 } .  We  will  also  comment  on  the 
linear  system  obtained  from  the  GFEM,  and  will  examine  Theorems  3.2  and  3.3 
in  the  context  of  a  classical  FEM  that  can  be  viewed  as  a  GFEM. 

Example  1:  The  classical  FEM  in  1-d,  based  on  continuous,  piecewise 
polynomials  of  degree  k,  is  same  as  the  a  suitably  chosen  GFEM.  We  show  this 
here  for  k  =  2  by  proving  that  the  finite  dimensional  approximating  space  used 
in  this  GFEM  is  same  as  the  classical  FEM  space. 

Suppose  Cl  =  I  =  (0,1),  and  for  a  fixed  positive  integer  N,  let  x3  =  jh , 
0  <  j  <  N,  with  h  =  1  /N,  be  uniformly  distributed  nodes  in  I.  We  consider 
the  “triangulation”  of  /  by  the  intervals  Ij  =  (xj,Xj+ 1).  The  standard  FEM 
space,  relative  to  this  triangulation,  is  given  by 

SFEM  =  {v  G  (7(0, 1)  :  v\j,  G  Vk(Ij),j  =  0, 1, . . . ,  N  -  1}.  (4.1) 

We  construct  a  GFEM  space  as  follows:  To  each  node  x3,  we  associate  a  function 
( t>j ,  which  is  the  usual  piecewise  linear  continuous  hat  functions  centered  at  x3 
such  that  <f>j{xi )  =  Sp.  We  let  u3  =  supp  <j>j  =  [xj-\,Xj+\\,  1  <  j  <  N  —  1.  For 
j  =  0,1V,  we  define  l6q  =  supp  4>o  =  [xo,a;i]  and  u>n  =  supp  4>n  =  [xn-i,  xn}. 
We  recall  that  the  sets  Wj' s  were  introduced  in  Section  2.  Clearly,  {<pj}jL0  form 
a  partition  of  unity  in  I  and  satisfy  (3.1)-(3.4).  For  the  local  approximation 
spaces  Vj,  0  <  j  <  N,  we  consider 

Vj  =  span{l,  x  —  x3}. 

We  then  define  the  GFEM  space  as 

N 

SGFEM  =  W  :  =  E  &(*)  (4-2) 

3=0 
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where 


h (x)  =  ai  +  fa ( x  -  xj)  e  vi -  a3 -  0i  G  R- 


The  functions  Ij  G  V3  are  only  defined  in  u>j,  but  since  <f>j(x)  =  0  at  Xj_ i  and 
Xj,  d>j(x)L(x)  has  a  natural  continuous  zero-extension  to  I.  We  will  show  that 

gFEM  __  gGFEM 


Since  the  functions  < f>j(x)lj(x )  are  continuous  on  /,  it  is  clear  that  functions 
in  SGFEM  are  also  continuous  on  I.  Also,  since  <pj  and  lj  are  piecewise  linear, 
it  is  clear  that  every  ip  G  gGFEM  a  £+),  piecewjse  quadratic  function.  Thus 
gGFEM  gFEM '  next  s}10w  that  gFEM  c  gGFEM ,  ie.,  for  a  given 

q(x)  G  gFEM ,  we  can  find  constants  a,;,/V  and  hence  U{x)  for  0  <i<N  such 
that 

N 

q(oc)  =  ^2  x  G  I.  (4.3) 

i= 0 


We  first  note  that  equality  of  q(x)  and  YiLo  <t>i(x)li(x)  the  nodes  Xj, 
0  <  j  <  AT,  implies 


N 

v(xj)  =  Mxj)h(xj)  =  o,-(.r;  )/,(./•.,  I  =  ay.  (4.4) 

i=0 

We  now  consider  the  function  YiLo  4>i(x)h(x)  with  these  cq’s.  Since  q(x)  and 
YiLo  Mx)k(x)  are  both  continuous,  have  same  values  at  nodes  Xj,  0  <  j  <  N, 
and  their  restrictions  to  the  IjS  are  quadratics,  they  will  be  equal  for  all  x  G  I 
if  they  are  equal  at  the  mid  points  of  Ij s,  i.e., 


N 

q(Xj+l/2)  =  X!  MXj+l/2)li(Xj+l/2),  0  <  j  <  N  -  1, 

i— 0 

where  Xj+1/2  =  Xj  +  h/2.  Imposing  these  conditions  yields 

N 


Q(Xj+ 1/2)  —  4>i(x j +  1/2)1  i(x 3  +  1/ 2) 


i— 0 


—  <l)j(Xj+l/2)lj  (Xj+1/2)  +  <t)j  +  l(Xj  +  l/2)lj  +  l(Xj  +  l/2) 

=  2  +  03(X3+l/2  ~  Xj)  +  aJ  +  l  +  /?7+l  (*J+l/2  -  *j+l)] 


a3  +  aj  +  l  +  2  (03  "  '^'+1  ) 


which  can  be  written  as 

&  -  0j+ 1  =  [^(v+i/s)  -  (atj  +  aj+i)]  0<j<N-l.  (4.5) 

For  an  arbitrarily  given  value  of  /?o,  we  can  solve  for  /V  1  <  i  <  /V  uniquely  in 
terms  of  /V  Using  these  /Vs  and  the  a,;’s  as  given  in  (4.4),  we  have  constructed 
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k(x),  0  <  i  <  N  such  that  (4.3)  is  satisfied.  Thus  SFEM  C  SGFEM ,  and  using 
the  fact  that  SGFEM  C  SFEM  (shown  above),  we  have  SGFEM  =  SFEM. 

It  is  well  known  that  for  k  =  2,  a  basis  of  SFEM  consists  of  nodal  hat 
functions  4>i(x),  0  <  i  <  N,  and  the  quadratic  bubble  functions,  Bpx),  0  <  i  < 
N  —  1,  given  by 


(x  ~  xi)(xi+ i~x),  Xi<x<xi+ n 


(4.6) 


0, 


otherwise. 


It  will  be  useful  later  in  this  section  to  have  an  expression  for  Bpx)  of  the  form 
(4.3).  From  (4.4)  with  q{x)  =  Bpx),  it  is  clear  that 


=  Bi(xj)  =0,  0  <  j  <  N. 


(4.7) 


Also,  since 


J  =  * 


Bi{Xj+ 1/2)  —  ^  4 

0,  j  ^  i, 


from  (4. 5), with  q(x)  =  Bi{x),  we  have 


solution  of  this  system  is 
Pj  = 


J  J.  i-i 

{  0. 

in  terms  of  /3o-  If  we  take  fi0  = 

l/h,  the 

p  1  <j<i, 
i  +  1  <  j  <  N. 

(4.8) 

we  get 

i 

'52<t>j(x)  (■ x~xj )• 

3= 0 

(4.9) 

The  above  expression  for  Bpx)  is  of  the  form  (4.3)  and  thus  Bpx)  is  a  linear 
combination  of  the  shape  functions  rjjk  of  SGFEM . 


Remark  f.l.  We  recall  from  Section  3  that  the  functions  in  the  local  ap¬ 
proximation  space  Vj,  for  j  for  which  ujj  fl  I’-|  ^  0,  must  satisfy  the  Dirichlet 
boundary  condition  on  ujj  f'\  T-| .  In  this  1-cl  setting,  if  the  exact  solution  u  of  a 
BVP  satisfies  the  boundary  condition  u(0)  =  0  at  x  =  0,  we  take  a0  =  0  and 

Vq  =  span  {2}; 
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the  functions  in  Vo  satisfy  the  boundary  condition  at  x  =  0.  Likewise,  if  it(  1)  =  0 
is  the  specified  boundary  condition  at  x  =  1,  we  take  a  at  =  0  and 

Vn  =  span  {x  —  1}; 

the  functions  in  Vn  satisfy  the  boundary  condition  at  x  =  1.  A  minor  modifi¬ 
cation  of  the  above  analysis  shows  that  SGFEM  =  SFEM  in  this  case  also. 

Example  2:  Consider  the  domain  Q  =  (0,  l)2  and  for  a  fixed  positive 
integer  N,  let  Xi  =  ih,  y-j  =  jh,  where  h  =  1/N  and  0  <  i,  j  <  N.  We  consider  a 
“triangulation”  of  by  the  squares  =  (xj,  x,+i)  x  ( yj ,  yj+i),  0  <  i,  j  <  N—l. 
The  nodes  of  this  triangulation  are  Ai}j  =  ( Xi,yj ),  0  <  i,j  <  N.  A  standard 
FEM  space  with  respect  to  this  triangulation  of  f l  is 

SFEM  =  {»£  C°(L 1)  :  u|Qi  .  G  Qk(nid)},  (4.10) 

where  Qk(flij)  =  span  {xlym}fm=0,  i.e.,  the  space  of  polynomials  of  degree  <  k 
in  each  variable.  It  is  possible  to  find  a  GFEM  space,  SGFEM ,  with  suitably 
chosen  partition  of  unity  functions  {<t>i,j{x,y)}  and  local  approximation  spaces 
Vij,  so  that  SGFEM  =  SFEM.  We  again  do  this  for  k  =  2.  For  k  =  2,  the 
functions  in  SFEM  are  C°  piecewise  biquadratics.  We  construct  a  GFEM  space 
as  follows:  To  each  node  we  associate  a  function 

=  (, h(x)<j>j(y ),  (4.11) 

where  <j>i(x)  and  are  one  dimensional  hat  functions  centered  at  Xj  and  yj 
respectively,  as  discussed  in  Example  1.  <j>ij  is  the  standard  piecewise  bilinear 
hat  function  centered  at  At  j  satisfying  =  1  and  c pij  is  zero  at  every 

other  node.  We  let  u>ij  =  supp  (j>ij  =  [xj_i,x»+i]  x  [yj_\,yj+\\.  We  note  that 
when  i  =  0  or  j  =  0,  we  replace  Xj_i  by  xt  or  yi-i  by  yt,  accordingly,  in  the 
definition  of  u>ij.  Similarly,  when  i  =  N  or  j  =  N,  we  replace  Xj+i  by  xt  or 
yi+i  by  2/i,  accordingly.  Then  {<pij}  satisfy  (3.1)-(3.4),  in  particular,  they  are 
a  partition  of  unity  on  O.  For  local  approximation  spaces  V)j,  0  <  i,j  <  N,  we 
take 

Vitj  =  span  {(x  -  Xi)l(y  -  yj)m ,  1  =  0,1,  m  =  0, 1}. 

Thus  Vij  is  the  space  of  all  bilinear  functions  defined  on  Ujj.  We  now  define 
the  GFEM  space  as 

N 

SGFEM  =  ^  .  ^(x,y)  =  ^2  <t>i,j(x,y)  h,j(x,y ),  (4.12) 

i,j= 0 

where 

li.  j (x.  y )  —  (i'ij ~\~ b'ij (x  Xi) -\- Cijiyy  ijj )  V d,j ( x  ) (y  yj ) ,  ,  c,j ,  d,j  G  M. 

We  note  that  lij  is  defined  only  on  u>ij,  but  since  4>i,j\dw  =  o,  , /,  ,  has  a 
natural  continuous  extension  to  f 1.  Thus  SGFEM  is  equivalently  given  by 
sgfem  =  span  {0.  ^  (j,  _  Xi)<j>itj,  (y  -  yj)(j>ij,  (. x  -  Xi)(y -  (4-13) 
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We  now  show  that  SFEM  =  SGFEM .  Since  the  functions  (j>i,jh.  j  are  contin¬ 
uous  in  £2,  it  is  clear  from  (4.12)  that  the  functions  in  SGFEM  are  continuous 
in  £2.  Also  since  lij  are  bilinear  on  each  rectangle  of  the  triangulation, 
the  functions  in  SGFEM  are  C°  piecewise  biquadratic  functions,  and  hence 
gGFEM  f-  gFEM  _  remajns  t0  show  that  SFEM  C  SGFEM . 

We  will  do  this  by  proving  that  every  element  of  a  basis  of  SFEM  is  contained 
in  SGFEM .  For  k  =  2,  a  well-known  basis  of  SFEM  consists  of  the  following 
functions,  which  can  be  grouped  into  four  categories: 

(a)  The  hat  functions  <j>itj(x,  y )  corresponding  to  the  nodes  A,;j,  0  <  i,  j  <  N. 

(b)  The  functions  S^j(x,y)  corresponding  to  the  line  segments  (Aij, 
0<i<N-l,0<j<N,  defined  by 

S-j  (x,y)  =  Bi (x)  ( y ) .  (4.14) 

Here,  Bi{x)  is  the  one  dimensional  quadratic  bubble  defined  in  (4.6).  We 
note  that,  for  1  <  j,  supp  =  [xi,Xi+i]  x  j/j  i..yy.ij]-  For  j  =  0,  the 
support  is  [xi,xi+ 1]  x  [yo,y\\. 

(c)  The  functions  S\  j(x,y)  corresponding  to  the  line  segments  (Aij,  AijJ+i), 
0<i<N,0<j<N  —  1,  defined  by 

V )  =  <t>%{x)Bj{y).  (4.15) 

We  note  that,  for  1  <  i,  supp  =  [xi-\,  Xi+i\  x  [j/j,  y,+i].  For  i  =  0, 
the  support  is  [x0,a;i]  x  [yi,yi+1\. 

(d)  The  functions  Biij{x,y ),  corresponding  to  the  rectangles  ,  0  <i,j  < 
N  —  1,  defined  by 

Bij(x,y)  =  Bi(x)Bj(y).  (4.16) 

We  note  that  supp  Bij  =  [xi,Xi+\]  x  [yj,yj+ 1]- 

It  is  immediate  from  (4.13)  that  c pij  £  gGFEM  £or  0  <  i,j  <  N.  Using 
(4.14),  (4.9),  and  (4.11),  we  have 

S\]]{x,y)  =  Bi(x)<!>j{y) 

1  1 

=  jl'52<t>i(x)(x-xi)My) 
n  1=0 

1  i 

=  Xi)<j>itj(x,y), 

1  1=0 

and  therefore  from  (4.13),  we  have 

S'jj  e  SGFEM,  for  0  <  i<  N  -  1  and  0  <  j  <  N. 
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Similarly,  using  (4.15),  (4.9),  (4.11),  and  (4.13),  we  have 

s-j{x,y)  =  *  J^(y  ~  yi)<t>iAx’y ) e  sGFEM, 
n  1=0 

for  0  <  i  <  N  and  0  <  j  <  N  —  1.  Finally,  from  (4.16),  (4.9),  (4.11),  and  (4.13), 
we  have 

Bij{x,y)  =  Bi(x)BJ(y) 

i  A ,  i  A , 

-^2<t>i{x)(x  -  xi)  -  <t>m(y)(y  -  ym) 

1= 0  J  L  m= 0 

i  j  r  x 

=  ^2^2  MA^miy)  J^{x  -  Xi){y  -  ym) 

l— 0  m= 0  *- 

i  j 

=  ( X  -  Xi)(y  -  yrn)4*l,m(Xi  y)  G  SGFEM , 

1=0  m—0 

for  0  <  *  <  N  —  1  and  0  <  j  <  N  —  1.  Thus  we  have  shown  that  all  the  basis 
elements  for  SFEM  belong  to  SGFEM .  Therefore,  SFEM  =  SGFEM . 

Remark  J^.2.  We  note  that  the  local  approximation  in  Example  2,  for 
the  indices  i,j  where  u>itj  flTi  ^  0,  can  be  chosen  such  that  all  lij(x,y )  € 
satisfy  the  Dirichlet  boundary  condition  on  Qij  D  r1?  i.e.,  lij(x1y)  =  0  for 
{x,y)  £  Co-ij  flT!,  and  Vii3  do  not  contain  constant  functions  for  these  indices 
i  and  j.  Moreover,  SGFEM  =  SFEM  for  any  k  in  (4.10),  and  thus,  in  this 
example  (also  in  Example  1),  the  GFEM  spaces  are  same  as  the  FEM  spaces 
corresponding  to  the  h-  as  well  as  p-  version  of  FEM.  We  further  note  that 
for  any  polygonal  domain  and  for  any  triangulation  of  Cl,  the  classical  FEM 
space  of  C°  piecewise  linear  polynomials,  can  be  viewed  as  a  GFEM  space  with 
standard  hat  functions  serving  as  the  partition  of  unity  functions,  and  where 
the  local  approximation  spaces  contain  only  constant  functions. 

Through  Examples  1  and  2,  we  have  shown  that  certain  classical  FEMs  can 
be  cast  in  the  framework  of  a  GFEM.  But  we  do  not  claim  that,  for  any  domain 
Cl ,  every  FEM  relative  to  every  triangulation  of  Cl  can  be  cast  in  this  framework. 
Our  main  reason  for  presenting  these  examples  is  to  illustrate  that  the  idea  of 
GFEM  is  a  generalization  of  the  idea  of  the  FEM. 

The  framework  of  a  GFEM  offers  more  freedom  in  choosing  shape  functions 
with  relatively  simpler  supports,  when  compared  to  classical  FEMs.  A  FEM  uses 
a  triangulation  of  the  domain  Cl,  or  a  mesh,  to  construct  piecewise  polynomial 
approximating  functions.  The  supports  of  the  shape  functions  (used  in  FEMs) 
are  union  of  “triangles”  relative  to  the  triangulation  or  the  mesh.  But  for 
domains  Cl  in  3-d,  with  complicated  geometry  ( e.g domains  with  voids  and 
cracks),  it  is  quite  difficult  to  generate  a  good  mesh  on  Cl.  One  of  the  important 
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aspects  of  the  GFEM  is  that  it  permits  the  use  of  partition  of  unity  functions  (in 
contrast  to  those  used  in  Examples  1  and  2),  whose  supports  may  not  depend 
on  any  mesh  ( e.g . ,  Shepard  functions  discussed  in  Section  2),  or  may  depend  on 
a  simple  mesh  that  does  not  conform  to  the  geometry  of  f 1  (see  [37]).  In  this 
sense,  the  GFEM  is  also  a  meshless  method  (see  [4])  and  this  feature  allows  us  to 
avoid  the  use  of  a  sophisticated  mesh  generator.  We  mention,  in  particular,  that 
for  the  partition  of  unity  functions  for  a  GFEM,  we  may  use  one  the  particle 
shape  functions,  e.g.,  RKP  shape  functions  (see  [26]),  used  in  meshless  methods. 
Another  important  aspect  of  GFEM  is  that  local  approximation  spaces  can  have 
functions  other  than  polynomials  (in  contrast  to  the  Vt  :j  used  in  Example  2), 
which  locally  approximate  the  unknown  solution  of  (2.1)  well.  Thus  the  shape 
functions  in  a  GFEM  need  not  be  piecewise  polynomials  (in  contrast  to  classical 
FEM),  and  the  approximating  functions  can  be  tailored  to  approximate  the 
unknown  solution  well. 

The  shape  functions  of  SGFEM  may  be  linearly  dependent  giving  rise  to  a 
singular  linear  system  (3.33).  This  can  be  easily  seen  in  Example  1  (k  =  2), 
where  there  are  2{N  +  1)  shape  functions  in  SGFEM ,  given  by  rjij  =  cf>i(x)(x  — 
Xi)i ,  j  =  0, 1,  0  <  i  <  N.  But  the  dimension  of  SFEM  in  (4.1)  with  k  =  2  is 
2N  +  1,  and  since  SFEM  =  SGFEM ,  we  have 

dim  SGFEM  =  dim  SFEM  =  2N  +  1  <  2{N  +  1). 

Thus  the  number  of  shape  functions  in  SGFEM  is  greater  than  its  dimensions; 
the  shape  functions  =  0,  l}.fl0  must  be  linearly  dependent.  Similar  con¬ 

clusion  is  also  true  for  the  shape  functions  of  SGFEM ,  given  by  (4.10),  in  Ex¬ 
ample  2  (also  see  [38]).  There  are  other  situations  in  which  the  shape  functions 
of  SGFEM  are  linearly  independent,  e.g.,  with  another  choice  of  partition  of 
unity  functions  as  shown  in  [28,  34],  But  the  shape  functions  could  be  “almost 
linearly  dependent”  giving  rise  to  a  severely  ill-conditioned  linear  system.  We 
will  discuss  the  solution  of  singular  or  ill-  conditioned  linear  system,  obtained 
from  GFEM,  in  Section  6. 

Finally  we  comment  on  Theorems  3.2  and  3.3  in  Section  3,  in  the  context  of 
the  FEM,  when  the  FEM  space  can  also  be  viewed  as  a  GFEM  space.  These 
theorems  are  fundamental  approximation  results  for  GFEM.  In  the  examples 
presented  in  this  section,  we  have  seen  that  SGFEM  =  SFEM ,  but  application 
of  these  theorems  on  SGFEM  does  not  yield  the  well  known  error  estimates  for 
the  FEM. 

In  Example  1,  the  FEM  approximating  space  SFEM  ((4.1)  with  k  =  2)  is 
the  space  of  C°  piecewise  quadratic  polynomials.  It  is  well  known  that 

||u  —  ufem\\s{Q)  <  Ch2\\u\\Hz(Q.),  (4-17) 

where  ufem  is  the  FEM  solution  relative  to  SFEM .  Here  u  is  the  smooth  (in 
H3(Q))  solution  of  an  elliptic  linear  Dirichlet  BVP  posed  on  fi  =  I  =  (0, 1)  with 
u(0)  =  w(l)  =  0.  Since  in  this  example,  SFEM  =  SGFEM ,  we  can  use  Theorem 
3.2  or  3.3  to  obtain  an  error  estimate.  Towards  this  end,  we  choose  £  Vv 
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0  <  j  <  N,  such  that 

ll«  -  $WeM  <  Ch\\u\\HHuj)  =  e2(j).  (4.18) 

Recall  that  Vj  =  span  {1,  (x  —  Xj)},  1  <  j  <  N  —  1,  Vo  =  span  {x},  and 
Vn  =  span  {(x  —  1)}.  Let  £“  =  <t>j(x)^‘j(x)  as  in  (3-11)-  It  is  easy  to  check 

that  (3.17)  and  (3.18)  hold  in  this  example,  and  thus  from  Theorem  3.3  and  the 
above  inequality,  we  get 

N 

llu-nl(n)  <  C5>2(j))2 

J=o 

N 

-  ch- \\u\\jj2(w.) 
i=o 

<  C/i2||u||^2(q). 

Thus,  using  (2.10)  with  ^  we  have 

\\u- uGfem\\£(q.)  <  ||w-|“||£(f2) 

<  Ch\\u\\H2(n),  (4.19) 

where,  ugfem  is  the  solution  of  (2.7)  with  S  =  SGFEM .  We  note  that  since 
gFEM  _  gGFEM ,  UpEM  =  ugfem-  But  (4.19),  which  the  based  on  Theorem 
3.3,  gives  only  O(h),  where  as  the  classical  estimate  (4.17)  gives  0(h2).  Thus 
Theorem  3.3  does  not  give  the  correct  order  of  convergence  in  this  situation. 
The  reason  for  this  loss  of  a  power  of  h  in  (4.19)  can  be  explained  as  follows:  The 
only  assumptions  on  partition  of  unity  functions  {4>j}  are  (3.1)-(3.4).  It  was  not 
assumed  that  {<f>j}  “reproduce”  linear  polynomials,  i.e.,  that  )CjLo  Xj<fij(x)  =  x, 
for  But  the  partition  of  unity  functions  {4>j}  used  in  Example  1  were  hat 

functions,  which  do  “reproduce”  the  linear  polynomials,  i.e.,  Y^J=o  xj(t>j(x)  =  xi 
for  x  £  I.  An  approximation  result  for  the  GFEM,  with  partition  of  unity 
functions  that  are  assumed  to  reproduce  linear  or  higher  degree  polynomials, 
will  be  reported  in  a  forthcoming  paper.  This  result  will  yield  an  0(h2)  error 
estimate  for  Example  1. 


5  Selection  of  Local  Approximation  Spaces 

As  we  have  seen  in  Sections  2  and  3,  the  local  approximation  spaces  play  a  cen¬ 
tral  role  in  the  GFEM.  We  discuss  the  selection  of  effective  local  approximation 
spaces  in  this  section. 

5.1  Selection  of  the  spaces  V]  using  the  available  informa¬ 
tion  on  the  solution  u 

As  mentioned  in  Section  3,  the  selection  of  local  approximation  spaces  Vj  is 
governed  by  the  available  information  on  the  exact  solution  u  of  Problem  (2.1). 
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In  this  subsection  we  discuss  some  types  of  available  information,  and  show  how 
it  can  be  used  in  the  process  of  selecting  Vj. 

(a)  The  available  information  on  u  is  in  terms  of  the  Sobolev  spaces: 

In  this  case  we  assume  that  the  only  available  information  on  u  is  that  it  lies 
in  Hm(ujj)  and 


(5.1) 


where  k  =  (. ki,k2),ki  >  0,  and  |fc|  =  k\  +  k2.  We  wish  to  select  the  spaces  Vj 
so  that 

sup  inf  ||u-£j||£(wn  is  small. 

u<^£  (ujj ) 

In  [3]  we  showed  that  if  we  know  only  (5.1),  then  the  space  of  polynomials 
of  degree  <  p  on  ujj  is  a  good  choice  for  Vj;  denote  this  space  by  Vj  =  W-p\ 
Then,  for  in  >  1, 

7  min(p,m—  1) 

E2(j)<CAf>  ,  (5.2) 

where  hj  =  diam  uij  and  C  is  independent  of  u,  h,p,  and  m. 

Remark  5.1.  The  estimate  (5.2)  is  the  best  possible  under  the  assumption 
that  the  only  available  information  is  (5.1). 

Remark  5.2.  From  (5.2)  and  Theorem  3.2  ,  specifically  (3.13),  we  obtain 
an  error  estimate  for  ugfem ■  Comparing  this  estimate  with  the  classical  FEM 
estimate,  we  see  that  we  loose  one  power  of  h.  This  is  because  we  have  assumed 
only  that  {4>j}  is  a  partition  of  unity,  i.e.,  that  it  reproduces  constants,  but 
possibly  not  linear  functions  (see  Section  4). 

(b)  The  available  information  on  u  is  in  terms  the  BVP: 

So  far  we  have  assumed  only  that  u  is  the  solution  of  the  BVP  (2.1),  i.e., 
we  know  nothing  other  than  that  it  satisfies  (5.1).  Often  we  know  more.  For 
example,  if  u  is  the  solution  of  (2.1)  with  a  =  1  and  /  =  0;  i.e.,  that 

(  A u  =  0,  for  (a:,  y)  £  O 

<  u  =  0  on  Tj  (5.3) 

l  ^=9onT2, 

then  u  is  a  harmonic  function.  Therefore,  in  this  situation,  we  use  harmonic 
polynomials,  instead  of  all  the  polynomials  in  WjP\  Let 

'HWjP'>  =  ju  e  :  v  is  harmonic  on  , 

the  left  superscript  H  denoting  harmonic.  Suppose  u):  is  star-shaped  with  re¬ 
spect  to  a  ball  and  duij  is  piecewise  analytic  with  internal  angles  ay  =  [IjH,  with 


23 


0  <  (3 j  <2  —  A,  A  >  0.  Then,  with  shape  functions  in  HWjP\  it  is  known  (see 
[29])  that 

,  ,  /losr»\ (2_A)(m_1) 

e2 (j)  <  CK^m)hm  1  ^  J  <j  p  >  m  -  1,  m  >  1,  (5.4) 

where  is  as  in  (5.1)  and  C  is  independent  of  u,  but  does  depend  on  the 

shape  of  u>j.  We  note  that  rate  of  convergence  in  p  depends  on  the  angle  of 
corners  of  the  boundary. 

Remark  5.3.  The  dimension  of  nW^p'>  is  2p  +  1,  whereas  the  dimension  of 

w!jp>  is  b>+1hjp+2)  _  Hence  for  a  given  asymptotic  rate  of  convergence,  the  space 
of  harmonic  polynomials  has  a  smaller  number  of  degrees  of  freedom  than  the 
space  of  standard  polynomials. 

Remark  5-4-  If  the  right-hand  side  is  not  zero,  then  we  have  to  add  additional 
shape  functions.  For  example,  if  /  =  1,  we  add  the  shape  function  £  =  x2  +  y2. 

Remark  5.5.  Because  there  is  a  known  relation  between  the  norm  ||tt||#m(w.) 
of  a  harmonic  function  and  its  trace  on  dujj,  we  can  express  (5.4)  in  terms  of 
an  appropriate  norm  of  u  on  duij . 

Remark  5. 6.  We  have  here  addressed  the  selection  of  shape  functions  for  the 
special  form  of  the  equation  in  (2.1),  namely  A u  =  0.  V.I.  Vekua  ([42])  and  I.N 
Bergman  ( [9] )  have  developed  a  theory  of  generalized  harmonic  polynomials  for 
differential  equations  with  analytic  coefficients,  i.e.,  functions  that  are  related 
to  the  differential  equation  as  are  harmonic  polynomials  related  to  Laplace’s 
equation.  For  a  discussion  of  generalized  harmonic  polynomials  in  connection 
with  the  equation 

A  u  +  k2u  =  0, 

see  [28]. 

Remark  5. 7.  Analogous  results  can  be  obtained  for  systems  of  PDEs,  e.g  the 
elasticity  equations,  and  higher  order  equations,  e.g.  the  biharmonic  equation. 

5.2  Selection  of  the  spaces  V)  when  u>j  has  a  complicated 
structure 

In  Section  5.1  we  tacitly  assumed  that  ojj  is  simply  connected.  Assume  now 
that  ujj  has  a  “circular”  hole  centered  at  some  point  (x,  y)  G  fl  and  consider  the 
problem  (5.3).  Suppose 

12  D  LUj  =  Uj ^  \  U>j2\ 

where 

uf]  =  {(x,y)  :\x-x\<h,  \y-y\  <  h} 

and 

=  {(A  y)  ■  (x  ~  A)2  +  (y  —  V )2  <  ^'2;  Where  5  <  h}, 
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and  assume  that  C  T2  and  g  =  0  on  dto^’ ,  ie.,  in  (5.3)  we  have  =  0 

( 2") 

on  duj\  .  We  consider  the  functions 


(2) 


Cj(r,  9)  =  (rl  +  r  l621)  sinl6,  l  =  1,  2, . . . 

;(2)  / 


£^(r,0)  =  (rJ  +  r  *<52/)  cos/0,  l  =  0, 1, . . . . 


(5.5) 


where  (r,  9)  are  polar  coordinates  with  respect  to  (x,y).  Clearly,  ,  i  =  1,2, 

are  harmonic  polynomials  satisfying  =  0  on  dui ^ .  Since  u  is  harmonic  in 
u>j,  it  can  be  expanded  in  an  infinite  (Laurent)  series  in  terms  of  the  functions 
in  (5.5): 


u(r,  9)  =  2«0  +  ]T  ai^(r,  9)  +  ^  9). 


(5.6) 


Thus  the  functions  in  (5.5)  can  be  used  as  shape  functions  and  linear  combina¬ 
tions  of  the  first  few  functions  in  (5.5)  provide  accurate  approximations  to  u. 


Because  ^  =  0  is  a  natural  boundary  condition,  which  need  not  be  explicitly 

(5.7) 


imposed,  we  can  use  the  functions 

rl  sin  19,  r~l  sin  19,  rl  cos  19,  r~l  cos  19. 


The  family  (5.7)  also  provides  accurate  approximations  to  u  on  luv 

Remark  5.8.  We  have  constructed  the  shape  functions  on  the  whole  plane 
with  one  hole.  If  the  domain  is  more  complex,  e.g.,  has  multiple  holes  as  in  a 
perforated  domain,  then  the  construction  of  the  shape  functions  is  more  compli¬ 
cated.  In  these  situations  we  can  use  (a)  numerical  construction;  (b)  analytical 
construction  based  on  conformal  mappings.  With  procedure  (b)  we  utilize  the 
facts  that 


1.  Conformal  mappings  preserve  the  harmonicity  of  the  functions;  and 

2.  Conformal  mappings  preserve  the  iL Cseminorm. 

Now  we  can  use  mapped  harmonic  polynomials  as  the  shape  functions.  For  a 
discussion  of  conformal  mappings,  we  refer  to  [22].  These  special  functions  are 
the  solutions  of  a  boundary  value  problem  on  the  domains  u>j  or  on  a  bigger 
domain  Qj  D  ujj .  We  call  these  problems  Handbook  Problems  because  they 
are  reminisant  of  the  handbook  problems  used  in  engineering.  These  problems 
(which  are  local)  can  be  solved  numerically  by  e.g.,  GFEM.  It  is  also  possible 
to  use  certain  analytic  formulas  similar  to  (5.6),  determining  numerically  the 
parameters  in  the  analytical  form  of  these  functions. 

So  far  we  have  assumed  that  u>j  is  a  domain,  i.e.,  a  connected  set.  In 
applications  the  GFEM  is  used  for  crack  propagation  problems.  Then  oj;/  is  “cut” 

by  a  line  into  two  domains  wj1'*  and  Uj  =  ujj1'1  U  and  fl  =  0. 
The  exact  solution  u  is  smooth  or  possibly  harmonic  separately  on  ct/-1 1  and 
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(2) 

u>j  ,  but  not  on  u)j  itself;  u  and  its  normal  derivative  are  discontinuous  across 

7  =  dLof]  n  dujf\ 

Here  we  have  to  create  the  space 


V,  =  vj->  = 


(1) 

on  lu  ~ 

(2) 

on  ujj 

so  that  there  is  a  =  (£  j1',  £;2^)  €  Vfp)  so  that 


Ap) 


\u-i 


■3  ’  s 3 

W||2 


fH £(wW)  18  SmalL 


(i) 

The  basic  Theorem  3.3  still  holds.  Denoting  by  x)  the  characteristic  function 
for  uJjl\  the  constant  function  mentioned  in  proof  of  Theorem  3.3  must  be 
replaced  by  {x^\xf^)-  Then  =  1,2  (respectively,  HW^  = 

=  1,2).  We  emphasize  that  in  vjp>  we  have  to  use  shape  functions 
in  ujj1'1  and  to ^  separately.  Then  we  get  analogous  results  as  before. 


5.3  Selection  of  the  spaces  Vj  when  the  solution  u  has  sin¬ 
gularities 

In  the  applications,  the  solution  of  (2.1)  can  be  singular  because  of  one  or  more 
of  the  following  reasons: 

1.  the  boundary  dfl  has  corners; 

2.  the  boundary  condition  changes,  e.g.,  from  Dirichlet  to  Neumann; 

3.  the  coefficient  a(x,y)  is  rough,  e.g,  it  is  piecewise  constant; 

4.  the  right-hand  side  is  not  smooth; 

5.  the  solution  has  a  boundary  layer. 

We  address  Items  1  and  2  only.  The  character  of  the  singular  behavior  of 
the  solution  of  (2.1)  is  well-known.  We  will  assume  that  the  boundary  dQ  has  a 
corner  at  A,  located  at  the  origin,  and  that  the  boundary  of  dQ  near  A  consists 
of  two  straight  lines;  this  assumption  is  only  for  the  sake  of  simplicity.  If  /  and 
g  in  (2.1)  are  sufficiently  smooth,  then  in  a  neighborhood  of  A, 

S 

u(r,  9)  =  ^  akrXk  logMfc  ripj(Q)  +  ( (r,  9),  (5.8) 

fc= o 


where  A^+i  >  A&,  Hk+ 1  >  Hk,  V’j(^)  is  a  smooth  function  of  9,  and  £(r,  9)  is 
smoother  than  any  of  the  terms  in  the  sum.  Here  (r,  9)  are  polar  coordinates 
with  origin  at  A.  We  note  that  (5.8)  is  also  true  when  Ti  (~l  T2  =  A,  which  is 
relavant  for  Item  2. 
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Now  we  select  the  shape  functions  in  Vj  to  be  the  functions  rAfc  logMfc  r*  (6*)  ^  k 
0, 1, . . . ,  s,  together  with  polynomials.  Then  the  error  of  the  approximation  of 
u  by  functions  in  Vj  is  only  the  error  in  the  approximation  of  £(r,  9)  by  polyno¬ 
mials. 

Remark  5.9.  There  is  a  large  literature  on  expansions  of  the  form  (5.8),  e.q.. 
[12,  20,  21,  31]. 

Remark  5.10.  An  expansion  similar  to  (5.8)  is  also  valid  for  elasticity  prob¬ 
lems. 

Remark  5.11.  If  we  have  g  =  0  in  (5.3),  then  Hk  =  0  in  (5.8). 

Construction  of  these  singular  functions  may  not  be  simple,  especially  in 
the  elasticity  problem.  Hence  a  numerical  treatment  is  unavoidable.  Either  we 
solve  the  associated  Handbook  Problem  (local)  problem  numerically  (with  the 
GFEM)  or  use  analytic  formulae  with  numerically  determined  parameters;  see, 
e.g.,  [33].  We  always  have  £  €  £(fi)  and  hence  it  is  not  necessary  to  use  the 
special  functions  in  (5.8)  as  shape  functions,  i.e.,  we  can  take  s  =  0  in  (5.8). 
However,  the  accuracy  when  using  only  polynomial  shape  functions  is  very  low. 

The  use  of  special  shape  functions  in  u>j  for  which  A  £  ulj  is  very  important. 
Also,  we  have  to  use  some  of  the  special  shape  functions  in  patches  uj:]  when 
A  ^  ujj,  but  ojj  is  close  to  A.  The  number  of  special  shape  function  needed 
depends  on  the  accuracy  requirement.  Determining  the  optimal  number  of 
terms  as  well  as  in  which  elements  special  shape  functions  are  needed  is  not 
simple.  Usually,  two  terms  in  patches  u)j  for  which  ulj  £  A  and  one  term  in  all 
u>j  that  are  the  direct  neighbors  of  these  patches  is  sufficient. 

5.4  Selection  of  the  spaces  V3  satisfying  the  Dirichlet  bound 
ary  condition 

If  ulj  flTi  =  0,  then  there  are  no  restrictions  on  the  approximation  functions  on 
dojj.  But,  if  ujj  fl  T-|  ^  0,  then  functions  in  Vj  must  equal  0  on  ulj  n  T^  Usually, 
it  is  not  difficult  to  create  such  functions.  For  example,  if  the  boundary  Ti  is 
a  straight  line,  or  a  circle  and  we  are  solving  the  Laplace’s  equation,  A u  =  0, 
then  it  is  easy  to  construct  such  functions. 

The  error  estimate  for  C2  then  depends,  as  before,  on  the  approximation 
properties  of  the  space  Vj . 

Remark  5.13.  If  the  Dirichlet  conditions  is  not  homogeneous,  functions  in 
Vj  must  satisfy  this  condition;  then  all  the  results  hold. 

Remark  5.14 ■  GFEM  constructs  uij  so  that  the  condition  \u)j  fl  T-|  > 

7  cliam  uij  is  satisfied.  This  is  easily  accomplished.  Then  there  are  no  diffi¬ 
culties  with  imposing  the  Dirichlet  boundary  conditions.  This  is  an  issue  with 
meshfree  method;  see  [4]  for  a  discussion  of  techniques  to  overcome  it. 
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6  Implementational  Issues  in  the  GFEM 

Implementation  of  the  GFEM  consists  of  four  major  parts,  namely: 

(a)  the  selection  of  local  approximating  functions; 

(b)  the  selection  of  partition  of  unity  (PU)  functions; 

(c)  the  construction  of  the  stiffness  matrix; 

(d)  the  solution  of  the  linear  system;  and, 

(e)  the  computation  of  data  of  interest. 

(a)  We  have  already  discussed  the  selection  of  local  approximating  functions, 
{£ji},  in  Section  5,  which  depends  on  the  available  information  on  the  unknown 
solution  u  of  the  problem  (2.1)  or  (2.6). 

(b)  The  primary  role  of  PU  functions,  {4>j},  in  GFEM  is  to  paste  together 

the  local  approximation  functions,  to  form  global  approximation  functions 

that  are  conforming,  i.e.,  global  approximation  functions  that  are  in  £n-  In 
theory,  any  partition  of  unity,  satisfying  (3.1)-(3.4),  will  suffice;  we  may  consider 
Shepard  functions  with  disks  as  their  supports,  as  described  in  Section  2,  or  finite 
element  hat  functions,  or  any  family  of  particle  shape  functions  used  in  meshless 
methods  (see  [4,  26]). 

But  the  choice  of  patches  {u>j}  and  the  associated  PU  functions  {<j)j}  affects 
many  aspects  of  the  implementation  of  GFEM,  e.g.,  (c)  and  (d).  We  first 
discuss  the  effect  of  patches  and  the  PU  functions  on  the  work  involved  in 

(c),  in  constructing  the  stiffness  matrix.  From  (3.34),  a  typical  element  of  the 
stiffness  matrix  is  of  the  form 

/  Vr]ik -VTjjidx.  (6.1) 

J  ojjC\u)i 

Since  these  integrals  are  evaluated  by  numerical  integration,  it  is  important  to 
choose  {u>j}  such  that  the  sets  {ujj  flwi},l  <  j,l  <  N  are  simple  domains, 
in  which  numerical  integration  could  be  performed  efficiently.  For  example,  if 
the  Wj’ s  are  disks  (in  M2)  or  balls  (in  R3),  a  typical  Uj  D  coi  is  a  “lens  shaped” 
domain,  and  accurate  numerical  integration  over  such  domains  is  known  to  be 
difficult.  We  note  however,  that  an  efficient  numerical  integration  scheme  for 
such  domains  was  reported  in  [13].  In  [34,  37],  u>j’ s  were  chosen  to  be  rectangles, 
and  a  typical  LUjClcui  was  also  a  rectangle.  It  is  much  easier  to  perform  numerical 
integration  on  rectangular  domains.  Thus  the  patches  {(Vj}  should  be  chosen 
so  that  the  sets  Wj  n  are  simple  enough  to  perform  numerical  integration. 
Moreover,  since  i)ri  =  the  integrand  in  (6.1)  has  terms  involving  {<fij} 

and  {V^j},  and  thus  the  numerical  evaluation  of  (6.1)  depends  also  on  the 
smoothness  of  the  PU  functions  and  their  derivatives  {V</>j}. 

The  choice  of  PU  functions  {<j)j}  also  affects  the  linear  system  (3.33).  We 
have  mentioned  in  Section  4  that  the  shape  functions  of  SGFEM  could  be  linearly 
dependent  or  independent,  depending  on  the  the  choice  of  PU  functions.  This, 
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in  turn,  leads  to  either  a  singular  or  a  non-singular  linear  system.  We  further 
note  that  the  condition  number  of  the  stiffness  matrix,  when  the  linear  system 
is  non-singular,  depends  on  the  choice  of  the  PU  functions.  Thus  the  choice  of 
PU  functions  affects  the  choice  of  the  linear  solver  used  in  (e),  since  the  choice 
of  linear  solvers  depends  on  linear  systems.  Finally,  the  constants  C\  and  Ci,  in 
(3.3)  and  (3.4)  respectively,  are  directly  related  to  the  choice  of  {<j)j},  and  these 
constants,  in  turn,  affect  the  constants  in  the  error  estimates  (3.13),  (3.14),  and 
(3.21).  We  note,  however,  that  it  may  not  be  wise  to  choose  the  PU  functions 
{<f>j }  based  only  on  any  one  of  these  effects.  The  choice  of  {4>j}  should  be 
balanced  with  respect  to  several  other  aspects  of  the  GFEM,  e.g.,  the  selection 
of  local  shape  functions. 

(c)  Evaluation  of  the  elements  of  the  stiffness  matrix  A,  in  (3.33),  involves 
more  than  just  ensuring  that  the  sets  {uj  fltc;}  are  simple  domains.  The  success 
of  GFEM  depends  on  evaluating  the  elements  of  A  with  high  accuracy.  Since  A  is 
symmetric,  only  the  upper  triangular  part  of  A  is  evaluated.  In  [37,  38],  the  same 
numerical  integration  was  used  simultaneously  to  evaluate  all  the  elements  in  the 
same  row  (the  diagonal  element  and  the  elements  to  the  right  of  the  diagonal 
in  the  same  row).  Also  numerical  integration,  based  on  adaptive  procedure, 
was  used  to  evaluate  these  elements.  In  the  problems  considered  in  [37,  38],  the 
diagonal  elements  of  A  were  always  dominant  and  a  low  tolerance  requirement  in 
the  adaptive  quadrature  for  evaluating  diagonal  elements  ensured  the  accuracy 
of  evaluation  of  off-diagonal  elements.  The  tolerance,  for  the  relative  error  in 
the  evaluation  of  the  diagonal  elements,  was  prescribed  as  0.01,  or  less,  of  the 
required  relative  accuracy  of  the  computed  solution. 

(d)  We  now  comment  on  solving  the  linear  system  (3.33).  We  have  mentioned 
before  that  the  stiffness  matrix  A  in  (3.33)  could  be  positive  semi-definite  or 
severely  ill-conditioned.  When  A  is  positive  semi-definite,  the  system  (3.33) 
has  non-unique  solutions.  We  have  mentioned  before  in  Section  3  that  the  lack 
of  unique  solvability  of  (3.33)  does  not  imply  that  the  GFEM  has  non-unique 
solutions. 

A  solution  of  (3.33)  can  be  obtained  with  (i)  a  specialized  direct  solver  based 
on  elimination,  or  (ii)  an  iterative  solver. 

(i)  The  linear  system  (3.33)  was  successfully  solved  in  [37]  using  the  direct 
method  of  multi-frontal  sparse  Gaussian  elimination  for  symmetric,  indefinite 
systems  that  was  developed  in  [17]  and  implemented  in  subroutines  MA47  and 
MA48  in  the  Hartwell  Subroutine  Library. 

(ii)  An  iterative  scheme  was  also  used  in  [37]  to  solve  (3.33),  which  we  de¬ 
scribe  here.  We  first  perturb  the  matrix  A  by  el,  where  e  >  0  is  small.  Let 

Ae  =  A  +  el. 

Clearly,  Ae  is  positive  definite.  We  first  compute 

c0  =  A~1b, 
r0  =  b  —  Acq  , 
zo  =  A^ro, 

Vo  =  Az0. 
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Then,  for  i  —  1,2,...,  we  compute 


1-  1 

Cj  =  cp+y^zj, 

3=0 

2- 1 

r*  = 

3=0 

=  Ae  Tii 

Vi  —  -T-2^2  ? 

until  the  ratio 

Icf^Cil 

is  sufficiently  small,  which  is  attained,  say,  for  i  =  I.  Then  cj  is  considered 
a  solution  of  (3.33).  In  practice,  we  have  seen  that  the  above  ratio  becomes 
sufficiently  small  in  one  or  two  steps.  For  a  numerical  example,  we  refer  to  [37]. 

(e)  Successful  solution  of  the  linear  system  (3.33)  yields  the  vector  c,  which 
is  used  to  compute  various  data  of  interest;  for  example,  approximation  of  the 
exact  solution  or  its  gradient  at  a  particular  point  x  €  fi.  This  data  is  obtained 
by  computing 


N  m(j)  N  m(j) 

ugfem(x)  =  EE  CjiVjiix)  and  S7uGfem{x)  =  EE  CjiV  r ]ji  (<r ) . 

3=1  *=  1  3  =  1  *= 1 

We  note  that  computation  of  r/j,;(x)  and  X7r)rl  (x)  involve  computation  of  <f>j{x), 
Zji(x),  V<j)j(x),  and  V£j,(a;).  There  are  other  data  of  interest,  e.g.,  stress  inten¬ 
sity  factors;  we  will  not  discuss  their  evaluation  in  this  paper. 


7  Applications,  Experience,  and  Potential  of  the 
GFEM 


We  have  discussed  the  basic  ideas  in  the  mathematical  foundation  of  the  GFEM 
in  the  simple  setting  of  linear  elliptic  BVPs. 

A  wide  variety  of  shape  functions  can  be  used  in  the  GFEM.  This  allows  the 
GFEM  to  successfully  approximate  non-smooth  solutions  of  BVPs  on  domains 
having  corners  or  multiple  cracks,  or  with  mixed  type  of  boundary  conditions 
-  Dirichlet  and  Neumann.  It  is  also  easy  to  construct  shape  functions  that 
are  smooth,  i.e.,  with  higher  regularity.  Thus  the  GFEM  can  be  used  to  solve 
higher  order  problems,  e.g.,  biharmonic  or  polyharmonic  problems.  Also,  the 
GFEM  with  smooth  shape  functions  can  be  used  in  problems  with  boundary 
conditions  involving  distributions,  in  which  situation  the  solution  of  the  BVP 
is  not  in  the  energy  space.  Furthermore,  the  capability  of  choosing  appropriate 
shape  functions  makes  the  GFEM  well-suited  for  solving  Helmholtz  problem 
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([28,  23])  and  certain  non-linear  problems  ([8]).  We  have  mentioned  before  that 
the  GFEM  either  does  not  employ  a  mesh  or  uses  a  mesh  only  minimally.  This 
allows  the  GFEM,  without  re-meshing  or  with  minimal  re-meshing,  to  be  used 
in  problems  involving  domains  with  changing  boundaries,  or  with  an  unknown 
boundary,  as  in  crack  propagation  problems  or  free-boundary  problems. 


Figure  1:  Example  of  a  perforated  domain 

The  GFEM  was  successfully  used  on  problems  with  complicated  domains  in 
[38,  40,  39]  using  simple  meshes,  and  thus  avoiding  complex  meshes  that  conform 
to  the  geometry  of  the  domain.  An  example  of  one  of  the  domains  considered 
in  these  papers  is  given  in  Figure  1.  We  note  that  the  voids  in  this  domain 
could  be  replaced  by  fibers.  In  fact,  the  positions  of  the  voids  in  Figure  1  are 
identical  to  the  positions  of  fibers  in  a  composite  material  and  were  obtained 
by  actual  measurement  ([2]).  Such  problems  were  successfully  solved  in  [39,  40] 
by  the  GFEM  using  simple  8x8  and  16  x  16  uniform  square  meshes  to  cover 
the  perforated  domain.  Detailed  analysis  of  the  accuracy  and  computational 
complexity  was  given  in  [40].  The  problem  with  perforated  domain  is  a  typical 
example  of  multi-scale  problems.  Moreover,  the  GFEM  was  used  on  problems 
with  boundary  layers  in  [16]. 

The  major  cost  of  the  GFEM,  when  applied  to  problems  with  complex  do¬ 
mains,  is  the  numerical  integration.  And,  as  mentioned  in  Section  6,  the  success 
of  the  GFEM  depends  on  efficient  numerical  integration  based  on  adaptive  pro¬ 
cedures.  Adaptive  numerical  integration  based  on  Simpson’s  rule  turned  out  to 
be  most  effective  in  the  problems  considered  in  [38,  40,  39]. 

The  ideas  in  the  GFEM  have  potential  of  being  used  in  other  frameworks. 
We  have  already  seen  in  Section  4  that  certain  FEM  approximation  spaces 
could  be  viewed  as  special  cases  of  GFEM  spaces.  Also,  the  approximation 
spaces  in  certain  meshless  methods  can  be  viewed  as  a  GFEM  space  (with 
constants  as  local  approximating  functions  and  the  particle  shape  functions  as 
PU  functions).  A  GFEM  space,  SGFEM ,  has  the  potential  of  being  used  in  the 
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context  of  mixed  formulations  of  elliptic  BVPs.  SGFEM  can  also  be  used  in  the 
framework  of  collocation  methods.  Of  course,  there  are  many  open  problems  of 
a  mathematical  nature  in  the  use  of  SGFEM  in  mixed,  collocation,  or  possibly 
other  methods.  The  problems  of  implementation  of  these  approaches  are  also 
open. 

The  effectiveness  of  the  performance  of  the  GFEM  (or  similar  methods)  on 
certain  benchmark  problems  has  been  shown  in  the  literature  [1,  25].  But  these 
benchmark  problems  are  so  simple  that  the  performance  of  the  classical  FEM 
on  these  problems  is  often  superior  to  the  GFEM.  The  future  of  the  GFEM 
or  other  similar  methods  is  uncertain  unless  their  superiority  is  established  on 
appropriate  realistic  benchmark  problems.  It  is  extremely  important  to  classify 
problems  where  these  methods  will  outperform  the  classical  FEM. 

Finally,  we  provide  a  list  of  problems,  where  the  GFEM  and  other  similar 
methods  have  great  promise  of  being  efficient  and  successful: 

•  Problems  with  non-smooth  solutions,  where  some  information  about  the 
solution  is  known,  or  could  be  obtained  by  a  local  numerical  computation. 
The  non-smoothness  of  the  solution  could  be  due  to  either  the  boundary, 
or  the  coefficients,  or  the  type  of  the  problem,  e.g.,  the  Helmholtz  problem. 

•  Problems  where  the  domain  is  so  complex  that  creating  a  mesh  by  a  mesh- 
generator  is  either  not  feasible  or  not  efficient.  We  note,  however,  that  a 
lot  of  progress  has  been  made  in  creating  efficient  mesh-generators  in  the 
last  decade. 

•  Problems  with  time  dependent  boundaries  or  free  boundaries  (i.e.,  prob¬ 
lems  with  unknown  boundaries).  Typical  examples  of  such  problems  are 
crack-propagation  problems,  seepage  problems  and  parachute  problems. 

•  Certain  non-linear  problems,  e.g.,  metal  forming  problems. 

8  Appendix:  The  Poincare  Inequalities 

In  this  appendix  we  outline  the  derivation  of  bounds  for  the  Poincare  constants 
C3  and  C4  introduced  in  Theorem  3.3.  These  bounds  will  be  in  terms  of  simple 
geometric  data  for  the  patches  u):r  Throughout  this  section,  x  and  y  denote 
points  in  M2. 

Theorem  8.1  Suppose  u>  is  convex,  d  is  the  diameter  of  u>,  and  10  contains  a 
disc  of  diameter  d  >  — .  Then 

f  /?\ 3/2  f 

IMIzdM  <  2«i  (  -  I  d  |M|£(w),  for  all  v  G  £(u>)  satisfying  /  avdx  =  0. 

(8.1) 


32 


Proof.  We  now  outline  the  proof  of  estimate  (8.1).  We  will  use  the  following 
result:  If  to  is  convex,  then 


||U  -  Va,s\\Ll(w)  < 


MH)1/2  fP\3/2 

\s\ 


d2\\v\\S(u),  for  all  v  £  H1^),  (8.2) 


where  S  is  any  measurable  set  in  to,  and 


va,s  = 


Ml 


[  avdx,  where  |S|a  =  f  ady. 
’la  Js  Js 


For  a  —  1,  this  result  is  proved  in  [19].  The  proof  of  (8.2)  is  a  mild  extension  of 
the  proof  of  (7.45)  in  [19].  Now  suppose  to  contains  a  disk  of  diameter  d  >  ^  . 
Then,  taking  S  =  to  in  (8.2)  we  get 


//3\ 3/2  f 

IMIl2(u;)  <  2/ci  I  -  I  d|MU(<«0>  for  all  v  satisfying  J  avdx  =  0, 

which  is  (8.1).  □ 

Let  to  be  an  open  set  in  R2  and  suppose  l  C  dto  is  an  arc.  For  x  €  w,  let 
Si(x)  =  the  convex  hull  of  {x}  U  l 
be  the  sector  subtending  l,  and  let  7 i(x)  be  the  angle  of  si(x). 


Theorem  8.2  Suppose  u>  is  convex,  d  =  diam  (a/),  and  Co  is  a  disk  of  diameter 
d  >  — ,  whose  closure  lies  in  to.  Suppose 


7 i(x)  >  70  >  0,  for  all  x  £  Co. 
(Such  an  oq  exists  since  the  closure  of  Co  lies  in  to.)  Then 


(8.3) 


\v\\li(u)  S 


<  <  (  - 


/  3/2  K27t1 

—  )  2/ci  +  —  )  -  /■  d  Hvllffa,'),  for  all  v  £  £(u>)  with  v\  1  =  0. 

\aj  \ a )  70  J 

(8.4) 


Proof.  We  now  outline  the  proof  of  estimate  (8.4),  which  is  in  two  steps.  We 
first  use  estimate  (8.2)  with  S  =  Co  to  get 


Ml  Ll(u)  < 


fi\  (71-lwl)1/2  2 

aJ  - - d  IMIU(^)  +  Ma,s||L2(W) 


<  I  - 

.  a 


1 3\ 3/2  (»H)1/2 


Next  we  estimate  ||v||  1.2(0,).  For  x  £  Co  and  y  £  l,  we  have 

r\x~v  I 
/o 


r\*-v  I 

v(x)  —  v(y)  =  —  /  Dr[v(x  +  r(cos0,sin0))]  dr, 
Jo 


33 


where 


y  =  y(9)  =  x  +  r(cos$,  sin#), 

(■ r,9 )  denoting  the  polar  coordinates  of  y  with  respect  to  x.  Now,  if  v(y)  =  0 
for  y  £  l,  we  have 

r\x-y\ 

v(x)  =  —  /  Dr[v(x+r (cos  9,  sin  9))]  dr. 

Jo 

Integrating  this  equality  with  respect  to  9  from  0  to  7 i(x),  we  get 

C7i0)  r\x-y(9)\ 


Hence 


v(x)j i(x)  =  —  /  /  Dr[v(x  +  r(cos#,sin#))]  drd9. 

Jo  Jo 

-7!(x)  Ax-vm  Dr[v(x  +  r(cos9,sin9))\ 


K®)l  = 


1 

li(x) 

1  r 


r  drdd 


<  — 


'0  Jo 

_  \Dv(v)\ 

70  Jsi(x)  \x-y\ 

1  f  \Dv(y)\ 


dy. 


7o  Ju  \x  -  y I 


dy 


=  -Vh(\Dv\)(x), 

7o  2 


(8.6) 


where 


(V^h)(x)=  f  |x-y|2(M  l)h{y)dy 
J  UJ 

is  the  Riesz  potential  of  h.  Squaring  (8.6)  and  integrating  over  u,  we  get 

IMU2(u>)  <  —  ll^i  (\Dv\)\\lH&)  <  —\\Vi(\Dv\)\\ L2(a/)-  (8.7) 

7o  2  7o  2 

We  have  the  following  estimate  for  the  Riesz  potential  from  Lemma  7.12  in  [19]: 


\V»h\\L2(u)  <  -TT1/2\uj\1/2\\h\\L2{oj). 
9 


Combining  (8.7)  and  (8.8)  yields 

IMUjco  < 


P\  7  nd 

—  F  £(«)• 
a)  70 


(8.8) 


(8.9) 


Combining  (8.5)  and  (8.9)  we  have 


IMUl(w)  < 


P\3/2  (ttIwI)1/2  2 11  11  f  M1/2  7TG? 

-)  pvj - d  M  £(u,)  +  (  ~  )  \~\i/2 —  u  £(«)•  (8.10) 

a  J  |w|  \a  J  \u\L/z  70 
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Finally,  since  w  is  a  disk  of  radius  d  >  AL,  we  get 

3/2 


MIl2K)  <  I  2«1  +  ^  j  d\\v\\s,  for  all  v  G  £(w)  with  v\i  =  0, 

(8.11) 


which  is  (8.4).  □ 

Remark  8.1.  In  Theorem  8.2  we  assumed  that  l  is  an  arc.  This  hypothesis 
can  be  considerably  weakened;  for  example,  it  can  be  a  disconnected  set. 
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